In [99]:
# Solving systems of equations

# To solve a system of equations, we will use the Python module
# 'sympy' and the function 'linsolve'.  We first have to tell Python which 
# variables we will be using. In this first trivial example, we will use "x" 
# as the variable.

# Import sympy:
from sympy import *
In [100]:
# Define our variable:
x = Symbol('x')
In [101]:
# Next, we define our equation(s).  It is assumed that the expression is 
# equal to zero.  For example, the RHS of the equation is x - 2 and the LHS
# is zero.
eqn1 = x - 2
In [102]:
# Now we can implement linsolve([...], [...]). In the first set of 
# square brackets, list the equations.  In the second set of square brackets,
# list the variables that you want to solve for.
soln = linsolve([eqn1], [x])
soln
Out[102]:
$\displaystyle \left\{\left( 2,\right)\right\}$
In [103]:
# Python returns the solution as a `FiniteSet'.  We can covert the FiniteSet
# to list via:
solnList = list(soln)[0]
solnList
Out[103]:
$\displaystyle \left( 2,\right)$
In [104]:
# Finally, we can index the list to get the final solution
solnFinal = solnList[0]
solnFinal
Out[104]:
$\displaystyle 2$
In [105]:
# Here's a less trivial system of two equations and two unknowns.
# First, we need to define a second variable.
y = Symbol('y')
In [106]:
# Here's the system of two equations and two unknowns.
eqn1 = 10*x - 3*y - 5
eqn2 = -2*x - 4*y - 7
In [107]:
# We now implement linsolve.
soln = list(linsolve([eqn1, eqn2], [x, y]))[0];
print('x =', soln[0], 'and y =', soln[1])

x = -1/46 and y = -40/23
In [108]:
# It is also possible to algebraically solve a system of equations in terms
# of other variables.  Below is a system of 6 complex equations and 6 unknowns.

# The example below solves for the currents in an all-pass filter.
# (Optional experiment in PHYS 231.)

# We must define the six currents that we wish to solve for as variables as
# well as the other relevant variables (R is resistance, C is capacitance,
# L is inductance, v is voltage amplitude, and w is angular freqeuncy.
v = Symbol('v')
R = Symbol('R')
w = Symbol('w')
L = Symbol('L')
C = Symbol('C')
In [109]:
# Here is a method to define all of six of the i0 to i5 symbols in a single line
i0, i1, i2, i3, i4, i5 = symbols('i0:6')
In [110]:
# Before we write our system of equations, first note that in Python you make
# a number complex by following it with a 'j' (no space).  For example,
# to enter z = x + j*y we would type 'z = x + y*1j'.  In these expressions
# j represents sqrt(-1).  Below we evaluate the square of of j which should
# be equal to -1.  
1j**2
Out[110]:
(-1+0j)
In [111]:
# To continue an input on a new line, use the backslash notation.
eqns = [i0 - i1 - i2, i1 - i3 - i4, i2 + i4 - i5,\
        i0*R + 1j*w*i1*L + i4*R + 1j*w*i5*L - v,\
        i0*R + i2/(1j*w*C) + 1j*w*i5*L - v, i3/(1j*w*C) - i4*R - 1j*w*i5*L]
In [112]:
# For some reason, 'linsolve' doesn't work here.  However, we can use 'solve'.
soln = solve(eqns, [i0, i1, i2, i3, i4, i5])
soln
Out[112]:
{i0: (I*C*L*v*w**2 + 2.0*C*R*v*w - I*v)/(2.0*I*C*L*R*w**2 + 2.0*C*R**2*w + 2.0*L*w - 2.0*I*R),
 i1: I*v/(-2.0*L*w + 2.0*I*R),
 i2: -C*v*w/(-2.0*C*R*w + 2.0*I),
 i3: C*v*w/(2.0*C*R*w - 2.0*I),
 i4: (-I*C*L*v*w**2 - I*v)/(2.0*I*C*L*R*w**2 + 2.0*C*R**2*w + 2.0*L*w - 2.0*I*R),
 i5: -v/(-2.0*I*L*w - 2.0*R)}
In [113]:
# As you can see the solutions is very complicated!  

# The output of 'solve' is an object called a 'Python dictionary'.
# The individual solutions can be accessed using soln[i0], soln[i1], ...
# Here is the solution for i0:
soln[i0]
Out[113]:
$\displaystyle \frac{i C L v w^{2} + 2.0 C R v w - i v}{2.0 i C L R w^{2} + 2.0 C R^{2} w + 2.0 L w - 2.0 i R}$
In [114]:
# We can use 'subs([],[],[],...)' to enter in numerical values for the various
# symbols.  I'll import 'numPy' too so that we can access pi. Here is a numerical
# value of i0.
import numpy as np
i0num = soln[i0].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)])
i0num
Out[114]:
$\displaystyle 3.75829862092635 \cdot 10^{-8} \left(2.51327412287183 + 0.579136704174297 i\right) \left(5026.54824574367 - 1158.27340834859 i\right)$
In [115]:
# We can force Python to tidy up the number using 'N()'.
N(i0num)
Out[115]:
$\displaystyle 0.0005$
In [116]:
# We can use an option with N() to specify how many sig figs to keep.  In the 
# output, the numerical value shows only 1 sig. fig. because it is exact.
N(i0num, 3)
Out[116]:
$\displaystyle 0.0005$
In [117]:
# Here are numerical values for all six of the currents.
print('i0 =', N(soln[i0].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i1 =', N(soln[i1].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i2 =', N(soln[i2].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i3 =', N(soln[i3].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i4 =', N(soln[i4].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i5 =', N(soln[i5].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))

i0 = 0.000500
i1 = 0.000194 - 0.000244*I
i2 = 0.000306 + 0.000244*I
i3 = 0.000306 + 0.000244*I
i4 = -0.000112 - 0.000487*I
i5 = 0.000194 - 0.000244*I
In [ ]: